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ABSTRACT 



Aims. Our objective is to study the vertical dust distribution in the circumbinary ring of the binary system GG Tau and to search for 
evidence of stratification, one of the first steps expected to occur during planet formation. 

Methods. We present a simultaneous analysis of four scattered light images spanning a range of wavelength from 800 nm to 3800 nm 
and compare them with (i) a parametric prescription for the vertical dust stratification, and (ii) with the results of SPH bi-fluid 
hydrodynamic calculations. 

Results. The parametric prescription and hydrodynamical calculations of stratification both reproduce the observed brightness profiles 
well. These models also provide a correct match for the observed star/ring integrated flux ratio. Another solution with a well-mixed, 
but "exotic", dust size distribution also matches the brightness profile ratios but fails to match the star/ring flux ratio. These results 
give support to the presence of vertical stratification of the dust in the ring of GG Tau and further predict the presence of a radial 
stratification also. 

Key words, radiative transfer - stars: circumstellar matter - methods : numerical - polarisation - scattering 



1. Introduction 

The first step of planet formation is the coagulation of dust 
grains, during low speed collisions (Dor mrnk et alj |2006, and 
references therein). Micrometer size grains grow to form aggre- 
gates that will finally give birth to kilometer size planetesimals 
and eventually planets. In parallel to grain growth, dust grains 
are expected to settle towards the disk midplane due to the con- 
jugate actions of the stellar gravity and gas drag. 

Dust settling arises because gas is pressure-supported (and 
thus has a subkeplerian velocity) whereas dust is not. Therefore, 
gas slows down the dust resulting in the dust settling towards 
the midplane and spiraling inwards to t he central star because 
of angular momentum cons ervation (see Weidenschil lindl 1 977t 
Stepi nski & Valageaslll996l) . This settling is highly dependent 
on the grain size with the larger grains (> 1 m) almost unaf- 
fected by gas drag and the smaller ones (< 1 /mi) so strongly 
coupled to the gas that they follow its motion. For the interme- 
diate si zes, the dust distribution dec ouples strongly from that of 
the gas (Barrie re^Fouchet et alj|2005h . 

The result is the formation of a dust sub-disk close 
to the equatorial p lane (e.g., ISafronov & Zviaginal [l969; 
Dubr ulle et al.ll 19951) . that may fragment and produce planetes- 
imals when the density of the dust layer exceeds the gas one 
dGoldreich & Ward 1973). Despite remaining questions and un- 
certainties - how dust grains overcome the "meter size bar- 
rier" without being accreted onto the star (e.g. Weidenschillin 
119771) or destroyed by high speed collisions dJones et al.ll 199' 
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Blu m et aT] |2000|) . how Kelvin-Helmholtz instabilities prevent 
the dust sublayer to fragment - models predict that dust grains 
can grow up to very large sizes and settle towards the midplane, 
leaving a population of small grains close to the surface and 
larger grains deeper in the disk. 

Such a stratified structure should have consequences on disk 
observables like their sp ectral energy distributions (SE Ds) and 
scattered light images dDulle mond & Dominikl [2004b . So far 
however, most studies of dust settling have been based on SED 
analysis only dApai et al.ll2004t iFurlan et al.ll2005l) . 

An alternative approach to probe the disk structure is through 
high resolution scattered light images. Such images of proto- 
planetary disks are becoming available over a wide wavelength 
range in a few cases, from the optical with HST, to mid-IR with 
ground-based adaptive optics. Analysis of these images with 
model fitting were usually restricted to a single wavelength but 
multi -wavelength studie s are becoming more and more common 
(see lWatson et al. 2006 and references therein). Because scatter- 
ing properties are highly dependent on grain size and because 
different wavelengths probe different depths into the disk, scat- 
tered light images over a wide range of wavelengths are poten- 
tially a powerful tool to detect dust settling. Indeed, attempts to 
simultaneously fit all wavelengths with a single model under the 
assumption of well mixed gas and dust indicate that the grain 
size distribution is likely more complex than that of the ISM 
(Mathis et al.l 1977b . In particular, a dependence of the maximum 
grain size with vertical distance above the disk midplane (ver- 
tical settling) sometimes appears to be needed dDuchene et al.l 
12001 hereafter D04). 
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GG Tau is of special interest to study dust set tling. Its geome- 
try, w ell known from millimeter emission maps dGuilloteau et al.l 
1999), is simple: the dust is distributed in a circumbinary ring 
between 180 and 260 AU and the inclination is known, around 
40°, allowing a wide sampling of the scattering phase function. 
Multiple wavelength studies of the disk dD04h have shown that 
I-band scattered light images probe sub-micron grains located 
about 50 AU above the disk midplane, whereas at longer wave- 
lengths, in the L' band, observations probe deeper regions of the 
disk (25 AU above the disk midplane) and reveal the presence 
of grains larger than 1 pm, strongly suggesting the presence of a 
stratified structure. 

In this paper, we investigate (i) whether multiple wavelength 
scattered light images of the GG Tau ring can be interpreted 
in terms of a disk model with a stratified structure and (ii) 
whether such a structure is required to explain the observations. 
In Sect. |2] we present model fitting of synthetic scattered light 
images. Two families of models are explored: (i) models where 
dust and gas are perfectly mixed in the disk and where optical 
properties of the grains are varied and (ii) models that include a 
stratified structure described in a parametric way. In Sect. [3] we 
explore if the model with a stratified structure, in better agree- 
ment with observations, is compatible with a more physical (hy- 
drodynamical) description of dust spatial distribution. The effect 
of realistic dust settling on scattered light images of the GG Tau 
ring is studied by coupling hydrodynamical simulations of the 
dust distribution with a radiative transfer code. We discuss our 
results in Sect. H] and conclusions are presented in Sect. [5] 



2. Parametric models of the ring 

2.1. Model description 

Model fitting is performed on the multiple wavelength scat- 
tered light images obtained in the Hubble Space Tele scope 
filters F814W, F 110W, F160\\Q by iMcCabe etail d2002l) and 
iKrist et all d2002h and on Keck- AO images in L' by ]D04 

Scattered light images are calculated by means of the M onte 
Carlo radiative transfer code MCFOST dPinte et al.l l2006). The 
general concept is to follow individual photon packets that ran- 
domly scatter in the disk, with probabilities given by the radia- 
tive transfer equations and dust optical properties, track them, 
and produce images when they exit the model boundaries. 

In this paper, we focus our analysis on the stratified structure 
of the ring and we adopt the same ring geometry and dust grain 
composition as lD04l The disk is assumed to be in hydrostatic 
equilibrium with a Gaussian vertical profile. The model parame- 
ters are the inner and outer rad ii (190 and 260 AU res pectively), 
the total dust mass (0.0013 M Q , Guillot eauet all 19991) . the verti- 



cal scale height (21 AU at 180 AU), the flaring index//(r) oc r 1U3 
and the surface density exponent £(r) oc r~ 17 . The density dis- 
tribution extends inside the inner radius with a radial Gaussian 
decrease with a 1 /2<? width of 2 AU : 



2(r) = E(190AU)exp 



1 [r- 190 AU\ 2 ' 



Am 



(1) 



with Aro = 2AU. The density is set to zero inside 180 AU. 

We assume dust grains to be homogeneous spheres (Mie 
theory) and adopt t he po rous dust grain optical properties from 
Mathis & Whiffen ( 1989, model A), with an 0.5 g.cm 3 average 



grain density. The grain size distribution n(a, r) is defined lo- 
cally, in each cell of the computation grid, and may then depend 
on the position within the disk, allowing to model disks with 
a stratified structure. Integrated over the whole disk, the grain 
size distribution is described by a power-law, dn(a) oc a p da. The 
minimum grain size is taken to be a mm = 0.03 pm and a max and 
p are considered as free parameters. With a max = 0.9 pm and 
p = -3.7, this dust model reproduces th e inter stellar extinction 
law up to 8 pm (Mat his & Whiffenlll989l) . lD04l showed that nei- 
ther a grain distrib ution typical of int erstellar medium nor the 
distribution used by IWood et al.l d2002l) to reproduce the SED of 
HH 30 are able to account simultaneously for both the scattered 
light images and I-band polarisation rates. In this section, we try 
to construct a grain size distribution that would reproduce both 
observed quantities simultaneously 

Models are monochromatic and are calculated at the central 
wavelength of the different filters (0.81, 1.0, 1.6 and 3.8 pm for 
the F814W, Fl 10W, F160W and L' filters). Each model is com- 
puted with a total number of 128 million photon packets. The 
pixel size is 0.04 "and images are convolved with a Gaussian 
kernel with a resolution similar to those of observations, i.e., 
FWHM= 0.065", 0.080", 0.128" and 0.091 "for I, J, H and 
L' bands, respectively. 

The quality of the models was determined by a x 2 min- 
imization method. The disk geometry wa s validated in pre- 
vious studies (e.g., iGuilloteau etai] 119991: ISilber et all 120001: 
iDuchene et alj|2004l) and in this paper we focus instead on the 
scattering properties, in particular the scattering phase function. 
In a first step, we only consider the azimuthal brightness profiles 
in our fitting procedure. These profiles were constructed with the 
same method as described in |D04| (see their section 4.2) using 
10° sectors. For each wavelength, we define a reduced^ 2 , based 
on the uncertainties in the images. Model uncertainties are neg- 
ligible with respect to the observations due to the large number 
of photon packets generated for each model. We finally define a 
total x 2 as the sum of the x 2 at each wavelength. From this total 
X 2 , a Bayesian analysis is performed to estimate the probability 
of occurrence of each parameter value. It also gives an estimate 
of the validity ra nge of e ach "best-fit" va lue (see for instance 
Press et aTl 1 19921 and also lLav et al.ll 19971 for an application to 
millimeter observations of HL Tau). Because our model is sim- 
plistic (the observations reveal systematic departures from the 
axysymmetry assumed in the model), this x 2 should not be in- 
terpreted in statistical terms but as defining a metric to compare 
the general behaviour of different modelfl 

2.2. Models without stratification 

Studies of grain growth processes show that grains with frac- 
tal dimension and/or with a size distribution that departs from a 
power-law, and whose scattering properties are probably notice- 
ably different from those of previously cited distributions, are 
produced. 

Modification of the grain size distribution or use of complex 
grain size distributions results in a dramatic increase of the pa- 
rameter space. Here, our aim is rather to verify how one can, by 
increasing in a reasonable manner the parameter space, modify 
the grain distribution used by D04 so that it reproduces simulta- 
neously all observations, without invoking any stratification. 



1 In the rest of the paper, we will refer to the F814W, Fl 10W, F160 
filters as I, J, H bands for commodity. 



2 The statistical interpretation of the Bayesian analysis remains 
however valid because systematic deviations between observations and 
models add a constant shift to the x 1 that does not modify the probabil- 
ity distributions. 
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Fig. 1. Azimuthal brightness profiles of the best model without 
stratification (solid line) superimposed to observed profiles (tri- 
angles). The model parameters are p — -5.5 and a max =100 fim. 
The position angle 0° corresponds to the disk semi-minor axis, 
on the front face of the ring. 



Scattered light images at multiple wavelengths probe differ- 
ent grains inside a single distribution. Thus, for a grain size dis- 
tribution typical of the ISM with a max = 0.9 fim, the median 
scatte rer ha s a size of 0.5 fim in the I band and 0.7 fim in the L' 
band dD04l) . 

A simple way to play on the characteristic size observed at 
each wavelength is to modify the relative fraction of each grain 
size via the slope of the size distribution. The goal is to introduce 
large grains, susceptible to reproduce the L' band image, and 
therefore larger than 1 fim, but in a small enough number so that 
their contribution at small wavelengths is negligible. 

To do so, we explore a 2 dimensional parameter space, by 
varying the maximum grain size between 1 fim (size slightly 
smaller than the size required to reproduce the L' band image 
in the work of |P04l) and 100/im (size large enough that the con- 
tribution from these grains is negligible on scattered light, even 
at 3.8 fim), and the index of the grain size distribution between 
-3.5 (as in the ISM) and -6.0, by steps of 0.5. 

Figure Q] presents azimuthal brightness profiles for the best 
model. Results are in agreement with observations at all wave- 
lengths with a global^- 2 of 57.28 (Table [TJ. The best model re- 
produces very well the images in the I, J and H bands but slightly 
overestimates the flux on the back side of the ring in the L' band. 

The Bayesian analysis shows that a single index for the grain 
size distribution, equal to -5.5, can reproduce the observations 



Table 1. x 2 of the best models without stratification. Each line 
gives the x 2 m the I, J, H and L' bands, as well as the total x 2 
defined as the sum of the individual ones for the models mini- 
mizing one of these x 2 , respectively. 



Model 


x 1 , 


xi 


xi 


2 

XI 


x 1 


Best model I band 


11.81 


16.70 


31.75 


25.77 


86.03 


Best model J band 


12.56 


16.35 


25.56 


3.47 


57.94 


Best model H band 


12.84 


19.23 


23.95 


21.45 


77.47 


Best model L' band 


19.66 


26.45 


55.62 


1.72 


103.45 


Best model 


12.16 


16.69 


25.99 


2.43 


57.28 
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Fig. 2. Bayesian probability distribution in the non-stratified 
case for grain size distribution slope (left) and maximum grain 
size (right). 




Fig. 3. Synthetic image of the GG Tau ring in the I band with 
only vertical stratification. The color map is in logarithmic scale 
and is inverted with bright regions in dark. The image size is 3.7 
arcsec. 

(Fig. [2j. While our crude sampling prevents us from accurately 
defining the allowed range for the slope index, we can readily 
exclude values smaller than -6.0 and larger than -5.0, leading to 
a tight constraint. Regarding the maximum grain size, we obtain 
that a minimum value of 3 fim is required but all values beyond 
that are equiprobable. 

2.3. Models with stratification 

We model the vertical stratification of the disk with a simple 
parametric law where the vertical scale height of the grains de- 
pends on the grain size: 



h Q (a) = h (a min ) 



(2) 



In the case of a perfect mixing between gas and dust, £ = and 
ho(a) - ho(a m i n ) is constant and corresponds to the gas scale 
height. 

The particular geometry of the GG Tau ring in which the 
back side of the inner rim of the ring is seen directly suggests 
that such a stratified model is too simplistic to describe scat- 
tered light observations. In Fig. [3] a synthetic map obtained in 
the I band with the stratification described above is shown. The 
inner edge of the ring, observed on the back face, presents a 
vertical brightness gradient (the surface is about twic e as bright 
as the mid plane) which is not seen in observations (Kris t et all 
20021120051) . 
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Fig. 4. Sketch of the vertical and radial stratification proposed to 
account for the scattered light images of GG Tau. 



This gradient is caused by the fact that, as a function of alti- 
tude at which photons hit the ring, they do not scatter on grains 
of the same size. Photons that scatter in the upper layers of the 
ring interact with small grains, which are isotropic scatterers and 
send a significant fraction of energy back to the observer. On the 
contrary, photons that scatter close to the midplane interact with 
larger grains, which are more forward-scatterers and send only a 
smaller fraction of energy back to the observeiQ. 

Thus, the different scattered light images cannot be explained 
only with vertical stratification. A complementary stratification, 
in the radial direction appears needed to explain the globally uni- 
form brightness observed on the back side of the ring. A config- 
uration where small grains would form a layer enclosing larger 
grains, as presented in Fig. |H might produce images in agree- 
ment with observations. At optical wavelengths, photons would 
scatter in the surface layer and meet small grains, whatever their 
height relative to the disk midplane. In the infrared, photons 
would penetrate deeper in the ring and scatter on larger grains. 

To represent such a stratification, we introduce an analytical 
description, similar to the one used for the vertical stratification, 
so that the density remains a continuous function. The charac- 
teristic length of the radial decrease of the surface density, Ar, 
becomes a function of the grain size, that we choose to describe 
with a power-law, whose index is the same as the vertical strati- 
fication one: 

Alia) = Ar ( — ) (3) 

\ ^min / 

with Ar fixed to 2 AU. 

We then explore a 2D parameter space, playing with max- 
imum grain size between 1 and 100//m as previously and the 
stratification index, that we sample between (no stratifica- 
tion) and 0.5. In this section, we set the index p of the integrated 
grain size distribution to a constant value of -3.5, following mod- 
els from lD04 

Figure [5] presents azimuthal brightness profiles obtained for 
the best model, compared to observed profiles. Results are in 
good agreement with observations at each wavelength, with an 
azimuthal profile similar at all wavelengths. The intensity on 
the back side of the ring, in the L' band is slightly larger, by 
about 1 cr, with respect to the observations. In order to quan- 
tify further the quality of the fit, we compare ;f 2 values obtained 
for the best model to those that we get by fitting one single 
wavelength (Table The global x 2 , 57.70, is comparable to 



3 Albedo effects are superimposed on those of the phase function. 
The albedo also varies vertically in the disk, but the relative differences 
are much smaller than for the phase function. 
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Fig. 5. Azimuthal brightness profiles of the best model with 
stratification (solid line) superimposed to observed profiles (tri- 
angles). The model parameters are £ = 0.15 and a max = 60 fim. 
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the one obtained from models without dust stratification (57.28, 
see Table [TJ. 

We note that, in the I, J and H bands, the best global model 
is only marginally poorer that the best models at each of these 
wavelengths but that the difference is more important for the L' 
band similarly to the models without stratification. 

Results of the Bayesian analysis are presented in Fig. [6] 
The stratification index is well constrained, with a probability 
curve peaking between 0. 125 and 0. 15 and a dispersion of about 
0.025. Regarding the maximum grain size, it appears that val- 
ues smaller than 5 fim are excluded but that beyond this limit, 
all values are roughly equiprobable, grains larger than 5 /mi con- 
tributing very little to the scattered flux at wavelengths smaller 
than 3.8 pm. 



Table 2. x 1 of the best models with stratification. The different 
lines have the same meaning as in Table Q] 



Model 


xi 


xi 


xi 


xi, 


Xlot 


Best model I band 


9.10 


16.67 


30.11 


2.44 


58.31 


Best model J band 


10.17 


14.95 


32.45 


7.93 


65.50 


Best model H band 


14.44 


20.54 


25.11 


32.13 


92.23 


Best model L' band 


17.09 


21.82 


48.20 


1.30 


88.41 


Best model 


10.18 


16.19 


29.03 


2.30 


57.70 
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Fig. 7. Observed (from Silb er et al]|2000l) and synthetic polari- 
sation at 1 pm as a function of the azimuthal angle. Polarisation 
maps have been binned by a factor 4, resulting in a pixel size 
of 0.172" and then divided in 20° sectors where the polarisation 
is measured, for each sector, in the pixel of maximum bright- 
ness. Black circles represent the polarisation extracted from ob- 
servations. The solid red line corresponds to the best model with 
stratification and a slope p = -3.5 and the dashed blue line to the 
model without stratification and a slope p = -5.5. The position 
angle 0° is defined as in the previous figures and corresponds to 
the disk semi-minor axis, on the front side of the disk. 



the variations of these ratios with wavelength. The model with 
dust stratification is in better agreement with observations: the 
flux ratio is roughly constant between the I and H bands and de- 
creases by about 15 % in the L' band. Even though it does not 
exactly reproduce the observations where the ratio slightly in- 
creases from the I to H bands, and then decreases in the L' band 
(=s 35 % smaller than in the I band), the global behaviour is re- 
spected, which is not the case for the model without stratification 
where a steady flux ratio decrease with wavelength is witnessed, 
with a factor around 6 between the I band and the L' band. 

Results presented in this section were based, for both cases, 
on the best model according to x 2 minimizations. The study of 
brightness profiles only gives a lower limit for the maximum 
grain size. Polarisation and colours show the same insensitiv- 
ity to the precise value of a max as the phase function. Results 
presented are therefore valid for all values a max > 5 pm and, in 
particular, cannot help to better constrain this parameter. 

3. Hydrodynamical modeling 

The parametric modeling of GG Tau's circumbinary ring pre- 
sented in the previous section points towards a stratified disk 
where both vertical and radial segregation of grains are needed. 

In this section, we aim at testing the validity of this hypoth- 
esis, using a more physical modeling of the processes of verti- 
cal settling and radial migration of dust grains. Results of hy- 
drodynamical simulations, where dust settling and migration are 
computed more realistically, are used as an input to the radiative 
transfer code to produce synthetic scattered light images. 



2.4. Polarisation and flux ratio 

Both models, with radial and vertical stratification or with a 
slope for the size distribution equal to -5.5, present equally good 
fits to the azimuthal brightness profiles. In order to discriminate 
between these two solutions, we compare their predictions of the 
polarisation and of the disk-to-star flux ratio with observations. 

Figure [71 prese nts the observed polarisation in the J band 
(Silb er etalJ boOO) as a function of the azimuthal angle com- 
pared with the polarisation predicted by the two best models and 
shows the differences between them. The model with stratifica- 
tion presents a polarisation rate around 15 % on the front side 
and peaking around 65-70 % on the back side of the disk whereas 
for the model without stratification but with a grain size distri- 
bution slope of -5.5, the polarisation rate is 10% on the front 
side and a maximum polarisation level of 80-90 % on the back 
side. Non e of the two mode ls provides a good fit to the obser- 
vations of Silb er et all d2000). The model with stratification per- 
forms slightly better (x 2 = 16 against ax 2 of 33 for the model 
without stratification) and reproduces reasonably well the polar- 
isation in the regions where it is the best defined, the front face 
of the disk. The model without stratification results in somewhat 
lower polarization rates (on the front side). However, both mod- 
els overestimate the polarisation on the back side of the disk. 
This behaviour is due to the high polarisability of the porous 
grains used here. An exploration of the dust properties, beyond 
the scope of the paper, is needed to obtain a better agreement 
with observations and, at this point, polarisation cannot be used 
as a strong constraint to discriminate between the two models. 

Table [3] gives all disk-to-star flux ratios for observations 
and models. These measurements depend on the ring geome- 
try, in particular to its projected surface on the plane of the sky. 
Because we did not explore this parameter, we focus here on 



3.1. Hydrodynamical code 

To compute the hydrodynamical evolution of gas and dust, we 
use Smoothed Particle Hydrodynamics ( SPH), a Lagrangian 
technique described by Monaghan] d 1992b . The equa t ions a nd 
approximations are rigorously established by Bicknell ( 1991). 

Our code was derived from Murray's code dMurravl [l 996) 
and is described in Barri ere-Fouchet et al.l (120051) . It is fully 3D, 
bi-fluid, non self-gravitating and locally isothermal with the tem- 
perature law given by T(r) cc r~ 3//4 , where T is the tempera- 
ture and r the distance to the star projected onto the midplane. 
The initial surface density is given by 2 oc r~ 3 ^ 2 , following es- 
timations of the Minim um Mass Solar Nebula (MMSN) and of 
lGuilloteauet"aT1 ( ll999l) . 

The dust phase, representing one single grain size at 
a time, is treated as a pressureless fluid and the coupling 
betw een gas and dust is done th r ough aerodynamic drag 
(see iMonaghan & Kocharvanl 119951: IMonaghan! 1 19971) . The 
drag exerted by dust on gas is neglected (see details in 
iBarriere-Fouchet et al.ll2005l) . 

In the initial state, dust and gas are considered as well mixed, 
which means that the dust disk is flared. This is correct for the 
smaller grain sizes when dust is strongly coupled to the gas, but 
for larger grain sizes, we can expect grains to start close to the 

Table 3. Disk-to-star flux ratios as a function of wavelength. 



Flux ratio Observations Models with Models without 



(%) (Duc hene et al. 2004) stratification stratification 

1 1.30 ± 0.03 092 057 

J 1.40 ±0.03 0.89 0.46 

H 1.54 ±0.03 0.91 0.30 

L' 0.97 ±0.09 0.77 0.10 
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Table 4. Simulation parameters for the circumbinary disk 



M, 2 



■■ 0.5 M 
: 0.65 M„ 



M eas = 0.13M o 
M dust = 0.01 M gas 
R iltik = 100 to 400 AU 
p d =1.0 gem -3 



a =35 AU e = 0.25 

* = 10, 7.5, 5, 2.5, 1 fim 

C s = C r- 3/s , C = 297.9 ms- 1 

S =S r- 3 / 2 

ttSPH = f = 0.1 

jSsph = 0.0 



midplane because they are produced by the growth of smaller 
grains that have already settled. Still, up to 1 mm, grains settle 
quickly enough to the midplane that this initial state does not 
really matter. The contribution of larger grains to the synthetic 
images is negligible (see Sect. |2.21 >. 



3.2. Parameters 

We run hydrodynamical simulations for a circumbinary disk 
around two stars of 0.65 and 0.5 M Q respectively, sepa- 
rated by 35 AU an d evolving on orbits with 0.25 eccentricity 
dDutrev et al.| [l994). We start from a disk with a 0.5 AU central 
hole and let the binary dig the hole during the 10000 first time- 
steps (=* 8 orbits at 100 AU) while the gas relaxes to a quasi 
equilibrium state before dust is injected. We then let the whole 
evolve for 40 000 more time-steps 32 orbits at 100 AU). 

The disk mass is 0.13 M with 1% dust in mass, its outer 
boundary is initially fixed to 350 AU and doesn't significantly 
spread beyond 400 AU. Each simulation includes about 25 000 
particles (gas+dust). The grain sizes range from 1 to 10 fim, we 
omit larger grains because their contribution to the scattered light 
images is negligible. Table [4] summarizes the simulation param- 
eters. 

The inner edge radius obtained in the simulations, about a 
hundred AU, is appreciably smaller than GG Tau's (180 AU). 
This is due to the uncertainty on the binary parameters for 
GG Tau and the depletion of the disk's central regions is 
therefore different. This simply introduces a scaling factor and 
does not alter in any way the results presented here. More re- 
cent binary parameters that could better account for the loca- 
tion of ring's inner edge have been obtained but still need to 
be confirmed by further ast rometric monitoring, as noted by 
iBeust & Dutrevl d2005Ll2006T) . 



3.3. Simulation results 

Figure [8] displays the dust and gas densities in the (r,z) plane 
at the end of the SPH simulations. The settling of dust appears 
clearly compared to the gas disk. The outer dust density contours 
show it to be more efficient for larger grains (10 /urn) than for 
smaller ones (1 fim) that remain strongly coupled to the gas. The 
radial distribution also differs from one grain size to the next. 
The 1 yum grains move closer to the star than the 5 and 10 /im 
grains (see also Fig. ITOb. which provides physical support to our 
parametric model of radial and vertical stratification (Sect.|2| in 
which the small grains form a layer surrounding the central parts 
of the ring where the large grains are present. 

The reason for this radial segregation is that larger grains re- 
spond more strongly to the pressure gradient and are more effi- 
ciently concentrated in the gas pressure maxima. Therefore, the 
radial extent of the density distribution decreases as the grain 
size increases (this is more easily seen in Fig. [TO}. Any asymme- 
try between the inner and outer pressure gradients will produce a 
shift in the peak of the density distribution according to the grain 




200 400 
r (AU) 



Fig. 8. Density contours for dust (left) and gas (right) at the end 
of the SPH simulations of the circumbinary disk for runs includ- 
ing 10, 5 and 1 /mi grains from top to bottom. 



size. In the case of GG Tau, the pressure gradient is stronger in 
the inner side than in the outer side, which results in the density 
peak being shifted outwards for larger grains. 

3.4. Interfacing the two codes: analytical fits to the data 

The SPH data are unfortunately too noisy to be used directly in 
MCFOST, even after smoothing on a regular grid. We therefore 
derive analytical fits that must consistently reproduce this differ- 
ential settling. To do this, we first perform a vertical fit of the 
dust density at the end of the SPH simulations at a given radius, 
the parameters of this fit are then fitted in the radial direction. 

3.4.1. Vertical fit 

We first decid e of a shape for the fit in the vertical direction. 
According to Garaud & Lin (2004), the time evolution of the 
dust density with height is self-similar, with initial conditions 
satisfying the hydrostatic equilibrium. This leads to the classical 
expression of density for a geometrically thin disk (in the verti- 
cally isothermal case): 



P(z) 



V2^H 



exp 



z 

27/2 



(4) 



Because of the self-similarity, the dust density profile remains 
Gaussian with a width and amplitude that evolve with time. 

We produce a vertical fit at a given radius ro with two param- 
eters p max and ctq so that 



P(Z) = Pmax exp 



•,2 \ 



2o-J 



(5) 



implying that H — <x and X = V27T ctq p n 



3.4.2. Radial fit 

For the radial fit, we use a parabola to fit the scaleheight H (see 
Fig. St): 



H(r) = Q2r +air + flo. 



(6) 
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Fig. 9. Fits for the 5 fim dust particles in the circumbinary disk, 
a. Scaleheight. b. Surface density. In both cases, the dotted line 
corresponds to the vertical fit only and the solid line to the sub- 
sequent radial fit. c. SPH density after smoothing on a regular 
grid. d. Vertical fit of the density, e. Complete fit of the density 
in both vertical and radial directions. 



We can thus reproduce a flared disk as well as a self shadowed 
one, as observed for the dust layer. This behavior of the aspect 
ratio has a strong infl uence on the SEP and the synth etic images 
as is highlighted by Dullem ond & Dominikl d2004l) and as we 
will see further on. 

In the case of a circumbinary disk, the surface density, plot- 
ted against the logarithm of the radius (In 2 = /(r)), looks like 
a Gaussian because the disk is actually a wide ring and there is 
a density drop towards the center as well as towards the outside 
of the disk. We choose a fitting formula where we look for the 
value of the three parameters E max , ro max , Cq' ver ify m g 



L exp 



' ( ln(r) - lnQom aJ) 
2< 



2 \ 



(7) 



The result of this fit is shown on Fig. [9j? for 5 fim grains. 

Finally, with the fits of H, of S and the theoretical expression 
of the density (Eq.|4]l, we get an analytical expression for it: 



P(r,z) = 



V2^ Zlo air 1 



x exp 



(ln(r) - Mrn max )) 2 
2 ^„ 2 



(8) 



This expression of the volumic density gives a much 
smoother distribution than the initial data (see Figs.[9j>e). 

The fitted scaleheights and surface densities for different 
grain sizes are plotted in Fig. [TUJ which clearly shows the dif- 
ferential settling, on the left panel, with the larger grains having 
the more impo rtant settling (because we are in the small parti- 
cle regime, see iBarriere-Fouchet e t al. (2005)). The right panel 
shows that the dust density maximum shifts towards larger radii 
as the grain size increases, and therefore as grains settle more 
efficiently in this size range (see Sect. 13. 3I >. 




100 300 300 400 500 50 

r (AU) 



200 300 400 500 



Fig. 10. Scaleheight (left) and surface density (right) for the 
small grains (1 to 10 fim) in the circumbinary disk. 



3.4.3. Interpolation between grain sizes 

Each hydrodynamical simulation only includes one grain size at 
a time. To synthesize scattered light images, we need to recon- 
struct a grain size distribution at each point of the disk. Because 
we do not consider interactions between dust grains and because 
all simulations, whatever the grain size used, predict consistent 
density for the gas phase, we can do this by superimposing the 
calculated densities for each grain size. The relative fraction of 
grains is fixed by assuming a grain size distribution, integrated 
over the disk, described as a power law with an index of -3.5. 
Finally, an interpolation, logarithmic with respect to grain sizes, 
is performed between the spatial distributions of grains resulting 
from the SPH calculations, to get a well sampled grain size dis- 
tribution. The smallest grains (a m ; n = 0.03 fim) are assumed to 
follow perfectly the gas distribution. 

3.5. Synthetic images 

Figure QT| compares the azimuthal brightness profiles in the I, J, 
H and L' bands resulting from our simulations of the circumbi- 
nary disk, with and without grain stratification. In the latter case, 
the aerodynamic drag of gas on dust is switched off and the dust 
distribution follows that of the gas for all grain sizes. The dif- 
ferences between models with and without stratification are es- 
pecially apparent on the back side of the disk. As for the ana- 
lytic stratification we used in Sect l2.3l we observe the smaller 
grains at the inner rim of the disk. They are isotropic scatterers 
and send back a larger fraction of light than the larger grains, 
of which we detect the scattered light when all grains are well 
mixed. We want to stress here that the models are not fits to the 
data, which are plotted as a reference in order to study the qual- 
itative behaviour of our models. 

The inclusion of dust segregation is a step in the right di- 
rection to explain GG Tau's observations. In particular, we note 
that the back side of the disk is particularly sensitive to the ra- 
dial stratification. The light coming from the front surface of the 
ring, which corresponds to angles between and 90° on one side 
and 270 and 360° on the other side, is in fact little affected by 
the segregation and both brightness profiles (with and without 
stratification) are almost identical. However, the back side is no- 
ticeably brighter in the stratified case, for the reasons mentioned 
above. 

It appears that the stratification resulting from the present 
hydrodynamical model is not sufficient to explain the observed 
brightness profiles. This may be due to the relatively short time 
on which the dynamical evolution is followed, or to physical pro- 
cesses which also lead to a spatial segregation of grains in the 



8 



C. Pinte et al.: On the stratified dust distribution of the GG Tau circumbinary ring 




Position angle (°) Position angle (°) 



Fig. 11. Azimuthal brightness profiles of the circumbinary disk 
with (solid line) and without (dashed line) stratification super- 
imposed to GG Tau's observed profiles (triangles). 



disk but are not considered here, like a preferred grain growth 
in the denser parts of the disk. Nonetheless, our hydrodynamical 
simulations demonstrate that it is physically plausible that the 
ring possesses an "inner wall" of small grains, as suggested by 
our parametric approach. 



4. Discussion 

The inclusion of stratification in a parametric model of GG Tau's 
disk (Sect. 12) allows to reproduce the scattered light images from 
the I to the L' band, in particular the azimuthal brightness dis- 
tribution which is remarkably similar at all wavelengths. This 
model is also in good agreement with the 1 yum polarisation 
data and globally reproduces the ring colours with respect to the 
central star. However, a simple parametric vertical stratification 
alone cannot explain all observations, especially the brightness 
of the inner rim. The presence of a radial segregation at the inner 
edge, with smaller grains closer to the stars than larger grains, 
also appears required to obtain synthetic images which repro- 
duce the observations. This kind of distribution can result from 
a preferred grain growth in the central, denser zones of disks 
where the collision probability is higher. It is also predicted by 
simulations of the dust dynamics around a binary system (see 
Sect. 13. 3b . which presents synthetic scattered light images show- 
ing a good qualitative match with the observations, as shown in 
SectEH 

Adding dust stratification and settling into models that al- 
ready contain several free parameters renders the problem of 
finding unique solutions even more acute. We explored a large 
parameter space and tried to quantify, through Bayesian analy- 
ses, the likelihood of occurrence of the models we presented. 
By doing so, we also found a solution with well-mixed dust 
that also matched the observed azimuthal profiles, at all wave- 
lengths. However, this solution appears physically less realistic 
than the expected one with settling because it requires a very 
steep slope for the dust size distribution. Yet, the solution ex- 
ists and we added observational constraints to confirm or reject 
its validity. When compared to polarisation data of the ring, the 
"well-mixed" solution proves only marginally compatible with 
the data, and certainly worse at matching observations than the 
model with stratification. More interesting, the well-mixed mod- 
els we explored are not able to reproduce the correct dependency 



of disk colours with respect to the star as a function of wave- 
length. 

The best "well-mixed" solution we identified is made of 
a distribution of spherical and homogeneous particles, with a 
slope of -5.5, noticeably steeper than th at of interstellar medium 
grains Models of dus t coagulation dDullemond & Dom inik 
120051: lOrmel et al.l l2007) show that, at least during the first steps 
of the process, the slope of the grain size distribution does not 
vary so much, in the size regime we are interested here. We 
succeeded until now in fitting a growing number of observable 
quantities by modifying the grain properties in a well-mixed 
model but this fit is still not perfect, and noticeably poorer than 
in the stratified case. These modifications are increasingly ar- 
bitrary and it becomes difficult to find physical causes to ex- 
plain them. Therefore, we consider this solution with a steep 
dust size distribution more "exotic" and less likely to occur. The 
model with stratification provides a better match to more data 
sets. Furthermore, stratification appears a more natural outcome 
of the disk evolution processes and for now it is our preferred 
mo del for the ring of GG Tau. 

iRrist et all d2005l) show evidence that the azimuthal bright- 
ness profiles change with time with front/back amplitude varia- 
tions around 20 %. These fluctuations remain moderate and tak- 
ing them into account only slightly affects the major observa- 
tional property of the ring in scattered light : the similar mor- 
phology of the disk over a large range of wavelengths. Therefore, 
our general conclusion that dust stratification is required to in- 
terpret scattered light images still holds but a more precise deter- 
mination of the degree of settling will probably require comtem- 
poraneous multiple wavelengths data sets. 

In all cases, the presence of grains larger than a few microns 
is required. Both models require a minimum value of a max of the 
order of 3 jum, indicating that the grains have evolved compared 
to those of the interstellar medium. The value inferred from our 
simulations assuming spherical grains is indeed a lower limit, 
a particle of complex shape being able to mimic the scatter- 
ing properties o f a smaller spherical and homogeneous particle 
(Rouleau 1996). This implies that grain growth towards sizes 
larger than one micron is happening in this circumbinary ring 
and that the first steps of planet formation are very probably be- 
ing initiated. 

5. Conclusions 

The definitive proof of a stratified structure in GG Tau's disk 
will probably rely on additional observations. However, strati- 
fied dust models are able to match several observable quantities 
of the GG Tau ring: in particular the azimuthal brightness pro- 
files from 800 nm to 3800 nm, but also polarisation measure- 
ments at 1 micron and the variation of star/disk flux ratios with 
wavelength. 

We have shown the attractiveness of polarimetric observa- 
tions to better constrain the grain properties. The broadening 
of the wavelength coverage also seems very promising, both in 
intensity and polarimetry. The detection of the disk at 10 /mi, 
where our models show that scattered light still dominates over 
the disk thermal emission, would enlarge the wavelength lever- 
age at our disposal and probe layers even closer to the disk mid- 
plane than the current 25 ALq The inclusion of the thermal 

4 Depending on grain models, and in particular on the fraction of sil- 
icates in the grains, the relative opacities at 3.8 and 10 yum vary widely 
and the zone that is probed at 10 fim is difficult to estimate before- 
hand. On the other hand, a geometrical study, similar to the one realized 
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emission allows additionally to derive an estimate of the ratio 
between the absorption and scattering cross-sections and to fur- 
ther restrict the families of solutions compatible with the obser- 
vations of a given disk. 

The modeling of the circumbinary disk in Sect. [3] allowed us 
to give a physical basis to the parametric modeling of GG Tau's 
ring in Sect. [2] Hydrodynamic models validate our hypothesis of 
a radial stratification of grains and are in qualitative agreement 
with the observations. It is indeed to the grains radial segrega- 
tion, rather than to their vertical stratification, that we are most 
sensitive in GG Tau's disk. The configuration we used does not 
exactly correspond to GG Tau's, due to uncertainties in the bi- 
nary's orbital parameters. The use of better parameters should 
allow us to make more quantitative comparisons with the data 
and to start to better constrain the processes causing the stratifi- 
cation in GG Tau's disk. 

We would like to stress here that validating the parametric 
approach by hydrodynamic computations is important because 
it allows the use of a faster tool to quickly estimate whether set- 
tling has occurred in a disk or not, before embarking on dedi- 
cated hydrodynamic calculations, which are more realistic but 
necessarily longer. 

On the long term, it will be necessary to include new physi- 
cal processes acting in addition to dust dynamics. Turbulence, in 
the calculations presented here, was only taken into account in 
a simplified manner. It would be interesting to see how a better 
treatment of turbulence limits the settling and affects the observ- 
able quantities. 

Grain growth is intimately linked to their settling. Because it 
enlarges the relative velocities, settling increases the number of 
collisions, which favours grain growth, which in turns increases 
the settling. The effects on the different observable quantities are 
potentially enhanced. The simultaneous modeling of both pro- 
cesses appears to be the next step towards a better understanding 
of the first stages of planet f ormation. It is cu rrently being imple- 
mented into the SPH code (lLaibe et alJ l2007) and the resulting 
simulations will be coupled to the radiative transfer code to con- 
tinue our study. 
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